function [rev, i1, i2, i3, rev_annual, mom, mom_skipping, cov12_23, cov12_34,  mom_annual, ref_s1, ref_s2, ref_s3, cov01_12, cov23_34, aaa, bbb] = main_func(A, v_theta, v_3, v_2, v_1, kappa, v_z, v_z3, mu, phi)
 
 % Just some preparation:
 
v_F3  = v_theta + v_3;
v_F2 = v_F3+v_2;
v_F1 = v_F2+v_1;

k_theta = (1/kappa)*v_theta;
k_F3 = k_theta+v_3;
k_F2 = k_F3+v_2;
k_F1 = k_F2+v_1;

%Date 3:

rho_3 = k_theta/k_F3;
k_theta_F3 = v_theta - k_theta^2/k_F3;
SIGMA_3 = k_theta_F3;

lambda_3 = A*SIGMA_3;
DELTA_3 = A^2*SIGMA_3;

%Date 2:

rho_2 = k_theta/k_F2;
k_F3_F2 = k_F3 - k_F3^2/k_F2;
SIGMA_2 = (k_theta/k_F3)^2 * k_F3_F2 + lambda_3^2/( v_z3^(-1) + DELTA_3);

 %SIGMA_3/( 1+ (A^2*SIGMA_3*v_z3)^(-1));

gamma_2 = lambda_3/( 1+ DELTA_3*v_z3);
lambda_2 = A*SIGMA_2;

%Date 1:

rho_1 = k_theta/k_F1;
k_F2_F1 = k_F2 - k_F2^2/k_F1;

SIGMA_1_ast = (k_theta/k_F2)^2*k_F2_F1 + (lambda_2 + gamma_2*mu)^2/ ( v_z^(-1) + A*( lambda_2+gamma_2 * mu^2)) ;
gamma_1 = lambda_2*mu + gamma_2*phi -  ( lambda_2+gamma_2*mu) * A* (lambda_2*mu + gamma_2*mu*phi)/ ( v_z^(-1) + A*( lambda_2+gamma_2 * mu^2));

LAM_1 = A* SIGMA_1_ast + gamma_1;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   theta,      F3,          F2       F1       z3,     z2,     z1,    	 
%   1           2           3       3       4       5       6     
vcov   = [...
    v_theta    v_theta 	v_theta v_theta 0      0     	0;...
    v_theta    v_F3         v_F3     v_F3     0      0		0;...
    v_theta    v_F3         v_F2     v_F2     0      0		0;...
    v_theta    v_F3         v_F2     v_F1     0      0 		0 ;...
    0          0           0       0       v_z3    0 		0 ;...
    0           0           0       0       0       v_z     0;...
    0           0           0       0       0       0       v_z];

%   theta,      	F3, 	F2       F1       z3,     z2,     z1,    	     
 P1 = [ 0 		0           0       rho_1   0 		    0  	                    LAM_1  ]';    
 P2 = [ 0       0           rho_2   0       0 	        (lambda_2+gamma_2*mu)   (lambda_2*mu+gamma_2*phi) ]';    
 P3 = [ 0    	rho_3     	0 	    0       lambda_3    lambda_3*mu 	        lambda_3*phi	]';    
 TH = [ 1 		0           0       0       0           0                       0]';

 aaa=rho_2;
 bbb=(lambda_2+gamma_2*mu);

 E1_P2 = [ 0    0   0  rho_2*(v_F2/v_F1)       0 	    0   (lambda_2*mu+gamma_2*phi) ]';    
 E1_P3 = [ 0    0   0  rho_3*(v_F3/v_F1) 	   0        0 	lambda_3*phi	]';    
 E1_TH = [ 0    0   0  1*v_theta/v_F1          0        0   0   ]';


R1 = P1;
R2 = P2 - P1;
R3 = P3 - P2;
R4 = TH - P3;

E1_R2 = E1_P2 - P1;
E1_R3 = E1_P3 - E1_P2;
E1_R4 = E1_TH - E1_P3;

ref_s1 = R1'*vcov*E1_R2;
ref_s2 = R1'*vcov*E1_R3;
ref_s3 = R1'*vcov*E1_R4;


rev = ( R1' * vcov * R2 +  R2' * vcov * R3 +  R3' * vcov * R4)/3;

mom  = (R1+R2)' * vcov * (R3 + R4);
mom_skipping = (R1+R2)' * vcov *  R4;

%cov02_23 = (R1+R2)' * vcov *  R3; 
%cov02_34 = (R1+R2)' * vcov *  R4;

cov01_12 = R1' * vcov *  R2;
cov23_34 = R3' * vcov *  R4;

cov12_23 = R2' * vcov *  R3;
cov12_34 = R2' * vcov *  R4;


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
SD1= (R1' * vcov * R1)^0.5;
SD2= (R2' * vcov * R2)^0.5;
SD3= (R3' * vcov * R3)^0.5;

sf = (2/3.1415927)^0.5;
Y_S = (SD1+SD2+SD3)/3 * sf;

SDM= ((R1+R2)'  * vcov * (R1+R2) )^0.5;
Y_L = SDM *sf;

scalefactor = 0.25*(1/2)^0.5 *8;

rev_annual =  rev/Y_S *scalefactor;

mom_annual   = mom/Y_L/2 *scalefactor;
%mom_skipping_annual = mom_skipping/Y_L *scalefactor;

i1 = R1' * vcov * R3; 
i2 = R2' * vcov * R4;
i3 = R1' * vcov * R4;




